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METHOD AND APPARATUS FOR SIGNAL RESTORATION 



This application includes a listing of a computer 
program. A portion of the disclosure of this patent 
document contains material which is subject to copyright 
protection. The copyright owner has no objection to the 
facsimile reproduction by anyone of the patent document or 
the patent disclosure, as it appears in the Patent and 
Trademark Office patent file or records, but otherwise 
reserves all copyright rights whatsoever. 

FIELD OP THE INVENTION 

This invention relates generally to the field of 
signal processing and particularly to the improvement of 
signal resolution, and the resolution of images represented 
by such signals, by compensating for distortions caused by 
the signal acquisition system. 

BACKGROUND OF THE INVENTION 

In the processing of various types of signals and 
images, including one dimensional signals (e.g., a function 
only of time) or multidimensional signals (e.g., an image 
signal which is a function of two or three dimensional 
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space) , it is often desirable to improve the quality of the 
signal from the signal acquisition system to correct for 
distortions introduced by the acquisition system. Such 
processing must be carried out in the inevitable presence 
of noise. Assuming a linear acquisition system, the output 
signal g from an acquisition system having an impulse 
response function h, responding to an actual signal f , and 
contaminated by random noise n, can be expressed as: 

g = h®f + n 

where ® represents the convolution operation. 

For many types of signal acquisition systems of 
great practical importance, the impulse response h is 
bandlimited in the frequency domain and varies over 
different portions of the signal, and cannot be calibrated 
accurately by empirical methods. Such acquisition systems 
include, by way of example, various types of microscopes, 
including widefield and confocal microscopes, scanning 
electron microscopes, scanning transmission electron 
microscopes, many other types of spectroscopic instruments 
including X-ray, infrared, etc., and one and two- 
dimensional gel electrophoresis processes. 

One area of particular interest is fluorescence 
microscopy, which permits the imaging of very weak light 
signals from luminescent molecules. In fluorescence 
microscopy, the detection optics of the microscope system 
is modified so that only certain wavelengths of light from 
excited molecules reach the detector, and the use of such 
microscopy with various fluorescent probes makes the 
technique very useful for imaging cells or parts of cells, 
indicating concentration of ions, in the staining of 
specific organelles, in detecting and quantifying membrane 
potentials, and in the detection of target molecules both 
inside and outside of the cell. 

In conventional fluorescence microscopy, the 
fluorescent intensity at a given point in a stained region 
can provide quantitative information about the target 
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molecules. However, it is well known that in any 
microscope when a point object is imaged in three 
dimensions, the point is imaged not only in the plane of 
focus but also in the planes around the focal plane. 
Furthermore, the intensity of the point falls off more 
slowly in the direction of the optical axis as compared to 
a direction perpendicular to the optical axis. Therefore, 
in a 3-D specimen being imaged, the intensity measured at 
any point in a given plane will be influenced by 
contributions of the intensities of neighboring points, and 
the contributions of points from out-of -f ocus planes have a 
greater effect on the quality of the image than 
contributions from in-focus points. In the spatial 
frequency domain, these distortions are manifested in two 
different ways. First, outside of a biconic region of 
frequencies in the three dimensional spatial frequency 
spectrum, the spectrum of the image is degraded by a strong 
low-pass function, and secondly, inside the biconic region 
of frequencies, all the frequency components of the image 
are removed during the image formation process. The 
severity of these distortions depends on the numerical 
aperture of the objective lens used for obtaining 3-D 
images. Higher numerical aperture lenses yield better 3-D 
images than lower numerical aperture lenses because of 
their superior optical sectioning capability. However, in 
addition to such optics-dependent distortions, image data 
is invariably corrupted by imaging noise. Thus, it is a 
significant problem to obtain good qualitative and 
quantitative information about the actual specimen from 
degraded images because if out of focus information is not 
excluded, unpredictable errors can affect both the visual 
and quantitative analysis of the specimen. 

A relatively recent development in microscopy is 
the confocal scanning microscope, in which only a single 
diffraction-limited spot in the specimen is illuminated at 
a time, with the light emitted from the illuminated spot 
being focused through an aperture which is "confocal" with 
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the spot of illumination- An image of a three-dimensional 
specimen is acquired by scanning a frame in the in-f ocus 
image plane and then moving the image plane to acquire 
additional frames of data until the entire specimen is 
scanned. In theory, in a confocal microscope the out-of- 
focus information can be removed, but in practice reducing 
the size of the confocal aperture results in a 
corresponding reduction in the signal-to-noise ratio. 

Thus, in both conventional and confocal 
microscopy, as well as in other signal processing systems, 
it is desirable to extract information from the original 
signal which would otherwise yield, when the signal is 
displayed as an image, a blurred, noisy image. The process 
of extracting unblurred images is called decon volution. 
Typical prior deconvolution processes require knowledge of 
the blurring function (impulse response) of the signal 
acquisition system (e.g., a conventional microscope or a 
confocal microscope) called the point spread function 
(PSF) . A technique utilizing projection onto convex sets 
and a form of Landweber iteration has been used for image 
restoration where the exact PSF is available through 
empirical measurements, e.g., by measuring the PSF of a 
microscope which involved imaging a diffraction limited 
point source. See Gopal B. Avinash, Computational Optical 
sectionin g Microscopy with Convex Projection Theory with 
Applications . Ph.D. thesis, 1992, University of Michigan, 
University Microfilms International, Inc., and Gopal B. 
Avinash, Wayne S. Quirk, and Alfred L. Nuttal, "Three- 
Dimensional Analysis of Contrast Filled Microvessel 
Diameters," Microvascular Research, Vol. 45, 1993, pp. 180- 
192. 

In most practical situations a complete knowledge 
of the PSF is not available, or the calibration procedures 
have to be carried out before each new specimen is imaged, 
which is often so time consuming as to make the procedure 
impractical. Thus, attempts have been made to develop 
procedures called blind deconvolution which do not require 
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complete prior knowledge of the PSF, and which attempt to 
simultaneously identify both the blurring function (the 
PSF) and the ideal unblurred image using the observed image 
data which are blurred and noisy. 

A special type of image formation occurs in 
microscopy where the image of a specimen is degraded by a 
band-limited PSF. In the paper by T.J. Holmes, "Blind 
Deconvolution of Quantum Limited Incoherent imagery: 
Maximum Likelihood Approach," J. Opt. Soc. Am. A, Vol. 9, 
1992, pp. 1052, seq, , a method based on maximum 
likelihood estimation for blind deconvolution was 
described. The success of the method was attributed in the 
paper to properly constraining the PSF estimate and the 
specimen estimate in the simulation studies. in the paper 
by V. Krishnamurthi, et al., "Blind Deconvolution of 2-D 
and 3-D Fluorescent Micrographs," Biomedical Image 
Processing and Three-Dimensional Microscopy, Proc. SPIE 
1660, 1992, pp. 95-104, blind deconvolution was applied to 
deblur three dimensional fluorescent micrographs from real 
specimens, a further form of this procedure was described 
in a paper by T.J. Holmes, et al., "Deconvolution of 3-D 
v:iu? Field and Confocal Fluorescence Microscopy Without 
Knowing the Point Spread Function," Proceedings of Scanning 
•93, 1993, III -57-59 , where the authors outlined the 
specific constraints on the PSF, on the basis that an 
empirical band-limit protocol has provided comparable 
results to deconvolution with the known PSF. These initial 
attempts to estimate images of specimens in the absence of 
complete knowledge of the PSF have several disadvantages. 
The method requires hundreds of iterations to converge to 
the final solution and, therefore, requires significant 
additional computation accelerating hardware to perform 
even a very modest size (64 x 64 x 64) 3-D image in a 
reasonable amount of time. The procedure constrains the 
PSFs using theoretically obtained bandlimit parameters, 
which generally are not appropriate. See, D. A. Agard, et 
al., "Fluorescence Microscopy in Three Dimensions," Methods 
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in Cell Biology, Academic Press, Inc., Vol. 30, 1989, pp. 
353-377. 

Blind deconvolution is also of interest in other 
applications of signal processing generally and image 
signal processing in particular. The application of 
maximum likelihood estimators using expectation 
maximization schemes for simultaneous blur and image 
identification in two dimensions has been actively pursued. 
See, e.g., M.I. Sezan, et al., "A Survey of Recent 
Developments in Digital Image Restoration," Optical 
Engineering, Vol. 29, 1990, pp. 393-404. Generally, these 
procedures have been computationally very intensive and 
require the use of computers with array processors to 
complete the processing in a reasonable amount of time. 
Other methods which have been proposed have the same 
disadvantages. See, e.g., B.C. McCallum, "Blind 
Deconvolution by Simulated Annealing," Optics 
Communications, Vol. 75, No. 2, 1990, pp. 101-105, which 
describes the use of simulated annealing, and B.L.K. Davey, 
et al., "Blind Deconvolution of Noisy Complex Valued 
Image," Optics Communications, Vol. 69, No. 5-6, 1989, pp. 
353-356, vhich describes the use of Weiner filtering for 
simultaneous blur and inage identification. 

SUMMARY OF THE INVENTION 

The present invention provides a rapid and 
accurate restoration of a noisy signal obtained from an 
acquisition system to simultaneously determine an estimate 
of the true or ideal signal and the response function of 
the acquisition system. Such acquisition systems can 
include optical imaging systems, such as conventional and 
confocal microscopes, which produce signals corresponding 
to an image in which the ideal signal would represent the 
actual unblurred image at a particular focal plane, and can 
be extended to three dimensional images of three 
dimensional specimens in which multiple planes of data are 
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utilized to produce a three dimensional array of data. The 
signal can also be a function of a four or higher 
dimensional variable and can be a function of a one- 
dimensional variable, e.g., a function of time or of linear 
direction. The present invention does not require prior 
knowledge of the response function of the acquisition 
system, so that there is no need to calibrate the 
acquisition system before acquiring a signal, e.g., before 
each specimen is measured in a microscope. The invention 
provides an extremely efficient process with convergence to 
a satisfactory estimate much more rapidly than previous 
processes, with fewer iterations required and fewer 
manipulations required for each iteration. 

The invention further allows the sensitivity of 
the signal processing system to be adjusted to the level of 
noise in the output signal of the acquisition system 
instrument. By estimating the noise variance in the output 
signal, the sensitivity of the signal processing of the 
invention can be adjusted to best suppress the effects of 
noise. Moreover, the invention can be applied to local 
overlapping segments of signal data under conditions in 
which the response function ic changing ovsr the signal, so 
as to obtain a locally accurate estimate of the response 
function and the ideal input signal for the particular 
segments of signal. 

If the acquisition system response function is 
changing over the entire set of acquired signal, the 
present invention first obtains the signal g from the 
acquisition system, stores the signal data in a memory, and 
utilizes a processing means in the signal processor to 
divide the signal data into smaller overlapping segments, 
and then processes each segment separately so that the 
response function h can be assumed constant for this 
segment of signal while over the entire acquired signal it 
may be variable. For each signal segment, the parameters 
needed for constraining the impulse response function and 
the signal function are first determined. If the response 
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function over the entire signal data set is substantially 
constant, the signal data need not be segmented, and the 
procedure can be carried out on the entire signal data set. 
The parameters include estimates of the background signal 
level, noise variance, and frequency bandlimits of 
frequencies of interest in the signal g. If the signal g 
is a function of a multidimensional variable, the 
parameters are determined with respect to each dimension 
and preferably along orthogonal directions. An iterative 
procedure is then started to estimate the ideal signal f 
using an initial selected estimate of the response 
function, preferably a unit impulse function. Constraints 
in the spatial domain are applied to the resulting estimate 
of f , which is then used to obtain an estimate of the 
response function h. Constraints are then applied to the 
resulting estimate of h in both the spatial domain and the 
frequency domain. The updated estimate of h is then used 
to estimate f , and the cycle is repeated until a selected 
criterion is met indicating that acceptable results are 
obtained. The estimates of f and h are updated in the 
frequency domain using the transforms of f and h, wherein 
the iteration drives the transforms of f and h toward 
convergence in the frequency domain. Where the entire 
acquired signal is divided into multiple segments of 
overlapping signal data, the iterative procedure is 
repeated for each of the individual segments, and the 
resulting estimates of the actual signal data for each of 
the segments are merged together to provide a global 
estimate of the actual signal. The apparatus of the 
invention may then display the estimated signal as an image 
to the user, for example, a 2-dimensional display on a 
video monitor of a microscope image, and, if desired, can 
display the estimated impulse response function, or the 
functions corresponding to individual segments where the 
original signal data are segmented. 

Further objects, features and advantages of the 
invention will be apparent from the following detailed 
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description when taken in conjunction with the accompanying 
drawings . 

BRIEF DESCRIPTION OF THE DRAWINGS 

In the drawings: 

Fig. 1 is a block diagram illustrating the 
present invention. 

Fig. 2 is a simplified perspective view of the 
components of a confocal microscope system in which the 
invention is embodied for exemplification. 

Fig. 3 is a flow chart of the major functions of 
the signal processing system of Fig. l. 

Fig. 4 A is a flow chart illustrating the major 
steps carried out by the signal processing system to 
estimate the impulse response function h of the acquisition 
system and the ideal signal f . 

Fig. 4B is a flow chart illustrating a variation 
of the sequence of steps of Fig. 4A which may be carried 
out by the signal processing system. 

Fig. 5 is a simplified diagram illustrating the 
spatial exteut of the result of deconvolution on two 
dimensional data. 

Fig. 6 is a simplified diagram illustrating the 
manner in which signal data can be segmented , processed by 
segment, and the results of the processing pieced together, 
illustrated for two dimensions. 

Figs. 7A-7F are graphs illustrating the 
application of the invention to a one-dimensional set of 
synthesized signal data. 

Figs. 8A-8D are graphs illustrating the effect of 
changes in signal-to-noise ratio (SNR) in the application 
of the invention to the signal data of Fig* 7B. 

Figs. 9A-9E are graphs illustrating the effect of 
changes in the selection of the bandlimit in the 
application of the invention to the signal data of Fig. 7B. 
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DETAILED DESCRIPTION OF THE INVENTION 

The present invention provides improvement of 
signal information which may be acquired from various 
sources, in the presence of noise, to reduce the effect of 
the noise and to compensate for the distortions introduced 
by the acquisition system. This process is illustrated in 
Fig- 1 in which a source 20 provides true information, 
corresponding to a true or ideal signal w f", to an 
instrument (acquisition system) 21 which has an (unknown) 
impulse response h. The acquisition system 21 also 
receives (or generates) noise n, and provides an output 
signal g which is a distorted version of the f corrupted by 
the noise n. A signal processing system 23 stores the 
signal g as digital signal data in a memory 24 and utilizes 
the signal data g to provide an estimate f of the actual 
source signal f and an estimate h of the acquisition system 
impulse response h. These estimates are provided to an 
output device 25, which may be an off-line storage device, 
computer memory, hard copy output, video display, flat 
panel display, and so forth. Typically, in accordance with 
the present invention, the signal estimate f is displayed 
on a display device as a two-dimensional image. 

There are many types of physical signal 
acquisition systems to which the present invention is 
applicable. Examples include optical microscopes including 
widefield and confocal microscopes, scanning electron 
microscopes, scanning transmission electron microscopes, 
many other spectroscopy instruments including X-ray, 
infrared, etc, and one and two dimensional gel 
electrophoresis. By way of example only, a specific 
application of the present invention is to a confocal 
microscope system shown generally in Fig. 2. The confocal 
microscope system includes an optical microscope 31 which 
examines a specimen placed on a stage 32. The specimen may 
comprise, e.g., biological cells, microcircuit components, 
etc. Confocal microscopy is especially advantageous in the 
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examination of biological structures such as cells which 
have three dimensional form. One of the advantages of the 
confocal microscope is the ability to focus at focal planes 
at various levels in the cell to provide a set of -slices- 
through the cell which constitute light reflected or 
transmitted from the focal plane at the selected level in 
the cell, m the microscope system of Fig. 2, the light 
transmitted through or reflected from the structures at the 
selected focal plane in the specimen constitute the source 
20 of Fig. i. This light is received by the microscope 31 
and transmitted via optical couplings 34 to a laser 
scanning confocal microscope analysis system 35. The laser 
scanning confocal microscope system 35 and the optical 
microscope 31 together constitute an acquisition system, 
which provides the output signal g in the form of an 
electrical output signal in digital form on a line 3 6 to a 
computer processor 38. The computer 38 can provide various 
types of signal processing on the signal data to provide 
displays of the signal data on a video monitor screen 40. 
Such laser scanning confocal microscope systems are well 
known. An exemplary confocal microscope system is the 
Odyssey* laser scanning confocal microscope sold 
commercially by Noran Instruments, Inc. of Middlston, 
Wisconsin. An exemplary suitable confocal microscope 
system is set forth in United States patent 4,8 63,226 to 
Houpt, et al , , confocal Laser Scanning Microscope, the 
disclosure of which is incorporated by reference. For the 
reasons described above, in such confocal microscope 
systems, the impulse response function may vary across each 
two-dimensional slice of image signal data through the 
specimen and also from slice to slice in the axial 
direction of the microscope. Further, some light 
originating from positions away from the pixel being 
scanned at any instant in the confocal microscope will be 
received by the detector in the microscope system and will 
corrupt the image signal data. The present invention 
compensates for such distortions, and also minimizes the 
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effect: of noise which may come from various sources, to 
provide two- and three-dimensional images for display on 
the video screen 40, or for use in other manners, in a much 
shorter period of time than has heretofore been possible. 

As noted above, the effect of the acquisition 
system response function h on the signal f from the source 
in the presence of noise n may be expressed as g « h®f+n, 
where ® represents convolution. The process to find the 
estimate f of the signal f is carried out in the signal 
processing system 23 of Fig. 1, which may comprise the 
computer 38 of the confocal microscope system of Fig. 2 
(e.g., a Silicon Graphics Indy workstation, an IBM AT-PC, 
etc.), a dedicated computer within the laser scanning 
confocal analysis system 35, or in dedicated hardware 
components. If the impulse response function h of the 
acquisition system is known to vary over the entire set of 
signal data for the output signal g, or is believed to 
vary, the signal data g can be segmented by the processor 
23 as indicated at the block 45 in Fig. 3 in a manner which 
is suitable for an overlap and save method. For confocal 
microscope signal processing, it is generally preferred 
thot if the foreground information is black, the image is 
inverted to represent the image in optical density units as 
indicated at 46 in Fig. 3. The signal processor then 
proceeds at 47 to estimate h and f given the signal g. 
When the estimation is completed, the processed signal 
segments are put back together and displayed as an image as 
indicated at 48. 

A more detailed flow diagram of the signal 
processing of the invention is shown in Fig. 4A. In 
carrying out the present invention to restore signals from 
bandlimited distortions, both spatial and spatial frequency 
domain constraints are used in iterative processing. As 
used herein, the term "spatial" refers to signals which are 
functions of one or more variables, wherein the variables 
may be space, time, a combination of both, or other 
quantities. In the processing of visual image data from 
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confocal microscopes, for example, the independent variable 
is a spatial vector in a three-dimensional array wherein 
each level in the array corresponds to one frame of data 
taken from a particular focal plane in the specimen. in 
the preferred processing of the invention, all of the 
parameters necessary for restoration are preferably 
computed from the observed signal data g. m the 
description of the processing set forth below, and shown in 
the flow diagram of Fig. 4A, lower case symbols represent 
signal entities (hereafter referred to as spatial domain 
entities) and upper case symbols represent the 
corresponding frequency domain entities. Transform 
operations are performed to transform spatial domain 
entities to the frequency domain and to inverse transform 
frequency domain entities to the spatial domain. Such 
transform operations can include Fourier transforms as well 
as other transforms, for example, Hartley and cosine 
transforms. The use of the Fourier transform is generally 
preferred because of well known characteristics and the 
availability of efficient computation processes such as the 
fast Fourier transform (FFT) . 

The process begins ds indicated at block 50 
wherein the input signal f„ is initialized to the observed 
data g, the number of iterations is selected (or a 
convergence criterion may be selected) and the impulse 
response function ho is initialized to a selected function, 
e.g., the Dirac delta function. An advantage of using the 
delta function is that its Fourier transform is a (known) 
simple uniform function of frequency, so that the initial 
step of transforming Ho can be eliminated, thereby saving 
one computational step. Moreover, the Dirac delta function 
satisfies a requirement of the PSF that the central pixel 
should have the maximum intensity. 

The signal data g is examined as indicated at 
blocks 51 and 52 to address the considerations discussed 
below. This examination may be carried out manually with 



WO 96/10794 



PCTAJS95/12259 



- 14 - 

the aid of computer processing or may be carried out 
automatically by the signal processor: 

(1) Spatial bounds on the signal. Where the 
impulse response function h is modeled as locally shift 
invariant but globally shift variant, processing is carried 
out on smaller segments of the signal. The signal 
sectioning may be carried out using well known overlap 
saving techniques. The edges of the signal are preferably 
padded to avoid inconsistencies that might arise when using 
circular convolution involving fast Fourier transforms 
(FFTs) . 

(2) Spatial frequency bounds on the impulse 
response function. These bounds are the most important 
constraints influencing the restoration of the impulse 
response function and hence that of the signal. Since it 
can be empirically shown that the theoretical impulse 
response function substantially differs from the 
empirically determined one, it does not make sense to use 
the theoretical parameters as has been done in the prior 
art. It is important to note that by restricting the 
parameters to theoretical values, one cannot restore a 
signal when spatially varying impulse response fuiict-.ions 
are involved. Instead, in accordance with the present 
invention, the spatial frequency bound for the unknown 
impulse response function is determined from the observed 
signal data g as described below. 

For a one-dimensional signal, the frequency 
spectrum of the signal g determined, for example, by 
performing an FFT on the signal data g, and the frequency 
spectrum is then plotted on a linear scale. For example, 
the plot of the frequency spectrum may be displayed to the 
user on the display device 25, which may be a video display 
terminal such as the video monitor screen 40 of Fig. 2. 
For two- and three-dimensional (2-D and 3-D) signals (or 
high dimensions) , the sections of the spectrum in mutual 
orthogonal spatial directions are determined and plotted. 
For example, for a 3-D signal the FFT may be taken for 
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signal data in x, y and z directions from a selected 
origin. In each of the plotted spectral sections, the zero 
frequency component (DC level or background) is preferably 
excluded. A characteristic pattern is generally evident 
from the plot, since it is the impulse response function of 
the instrument 21 and not the true signal which determines 
the bandlimit of the observed data. The extreme high 
frequency portion of the plotted spectrum is assumed to 
represent the underlying noise processes, and the cutoff 
frequencies in a signal in a given axis (e.g., along x, y 
and z axes for a 3-D signal) is determined by locating the 
point where the signal frequency content drops off to the 
noise content as one traverses on the spectrum from the 
high frequency end to the low frequency end. it should be 
noted that traversing from the low frequency end to the 
high frequency end using the above criterion can result in 
incorrect cutoff points because the spectral magnitude can 
have zeros. if the impulse response function is spatially 
varying, different bandlimits can be computed by sectioning 
the signal data as discussed above. Although there is no 
generally applicable theory which can be used for 
dstsirmini,-!Q proper sectioning of a signal for which th*» 
response function varies, the signal can be sectioned into 
relatively small segments of a convenient selected size and 
then processed by determining individual bandlimits on the 
assumption that the response function is invariant over 
each small segment, since the impulse response functions 
typically vary slowly across the entire signal, this 
assumption is usually valid. 

(3) Spatial extent of the impulse response 
function. The spatial extent of the unknown impulse 
response function is also unknown. In the present 
invention, an upper bound on the spatial extent of the 
impulse response function is used. This upper bound on the 
spatial extent is preferably usually overestimated by 
empirical determination of blur extent around edges or 
sharp transitions in the signal in each dimension. Since 
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the spatial extent of the truncated, bandlimited impulse 
response function is dependent upon the spatial frequency 
bound, one can use the spatial frequency bound to control 
the spatial extent as well. In cases where spatial extent 
is deliberately set to be larger than the one required, the 
impulse response function values are zero or very close to 
it provided the frequency bound is chosen properly as 
described above* 

(4) Noise variance bound on the signal. The 
iteration procedure requires a parameter, designated herein 
as X, which is computed from the noise variance bound on 
the signal. An upper limit for the noise variance can be 
computed either interactively or automatically. For 1-D 
signals, it is easier to interactively compute the noise 
variance by selecting a segment of the signal g where the 
variations are known to be due to noise only. A similar 
approach can be taken to compute the noise variance bound 
for 2-D signals. However, for well sampled 3-D signals, in 
which there are slight variations between the adjacent 2-D 
sections, an upper bound can be automatically determined by 
computing the variance of the difference sections. 
Difference sections are obtained by subtracting the 
constituting 2-D sections pairwise. 

Estimation of the ideal signal and the impulse 
response function is then performed iteratively using a 
selected form of Landweber iteration for both the signal 
and the impulse response. The spatial frequency 
bandlimits, discussed above, are used to constrain the 
frequency domain estimate of the response function during 
each iteration, and constraining the estimates in this 
manner, with proper choice of the bandlimits, facilitates 
the processing of signals in accordance with the present 
invention much more rapidly than prior processing systems. 
This iteration is preferably and conveniently carried out 
in the signal processor 23. 

Each iteration follows the following sequence (at 

iteration number k) : 
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(1) given the existing estimate in the frequency 
domain of f , F k , and the existing estimate in the frequency 
domain of h, H*, determine an improved estimator T k as 

T k » F k + DH/ (G-H k F k ) (1) 

where G is the frequency domain representation of the 
5 signal g and D is a bounded linear operator chosen such 

that ||DH k H k ||<2, and W * M indicates complex conjugate. 
Several choices for the operator D are available. A 
preferred operator is such that 

dh; = (2) 

l+XH k H k 

where X is the selected parameter discussed above. It is 
10 seen that as the estimators and F k approach the true 

functions, H k F k should approach G, and (G - HfcF k ) should 
approach zero. 

(2) F k is transformed to the spatial domain and 
the result is constrained to have real and positive values 

15 inside the selected spatial extent and zero outside to 

yield the updated estimate of f , f k * r . 

(3) The borders of f k+ , are modified a.«; 
conventional to take into account the periodic properties 
of the discrete Fourier transform. 

20 (4) f k+l is transformed to the frequency domain 

to yield the updated estimate in the frequency domain F k+l . 

(5) The frequency domain estimate of h is then 
updated as an improved estimator H k in accordance with: 

H^H^DF^CG-H^,) (3) 

where D is again a bounded linear operator. A preferred 
25 operator is 
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DF ^»_±f£_ (4) 
l+XF kM F kM 

It is seen that as the estimates improve, (G-H^F^,) should 

approach zero. 

(6) The chosen frequency domain bandlimits are 
applied to H k to eliminate frequency components above the 
selected bandlimit in each spatial frequency axis, the 
result is transformed to the spatial domain, and the result 
in the spatial domain is constrained to be real and 
positive inside the selected spatial extent and zero 
outside, and then the data values are normalized — for 
example, so that the sum of all values constituting the 
updated estimate h k +i is equal to one. 

(7) If the iterations have not been completed, 
the updated estimate h^ is transformed to the frequency 
domain to provide H^, and the iterative procedure above is 
repeated* 

The iterations are completed and the final 
estimates of f and h are made when a selected criterion 
reached. For example, this criterion may be simply that a 
prechcsen number of iterations have been complete*, or it 
may be related to the amount of change of the estimates 
during an iteration. 

An implementation of this iterative processing is 
illustrated in Fig. 4A. The signal update is performed at 
block 53 using the relation: 

XH; ( G-H k F t) > (5) 
k l+XH k # H k 

where k is the iteration number and H* is the complex 
conjugate of H. As discussed above, the parameter X may be 
selected based on empirical considerations from examination 
of the signal g or from prior experience with similar data 
or instruments, or it may be computed as the Lagrange 
multiplier of the minimization of the following convex 
constraint during the first iteration: 
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minimize {|r-f p || 2 , 
subject to Hg-(h0T p )D 2 = * v (6) 

N-l 

where for an N-vector x, |xfi 2 * £x*, d v is the upper bound 
of the estimated noise variance "and f p is the projection of 
the initial estimate of the signal onto a set containing r 
such that 

Dg-(h®r> || 2 £S V (7) 

The same X is used throughout the procedure, and the X is 
preferably computed in such a way that the noise variance 
is overestimated to avoid any problems with the restored 
signal stability. 

The estimate in (5) is then updated (block 54) 
through the sequence: 

* Take the inverse Fourier transform of F k and 
constrain the result to be real and positive inside the 
selected spatial extent and zero outside to yield f k+1 

* Then, if the selected number M n w of iterations 
has not been completed (block 55) , modify the borders of 
f k+ , to take into account the periodic properties of the 
discrete Fourier transform and take the Fourier transform 
of f k+1 to get F k+1 (bloc* 5C) . 

The impulse response function update is then 
carried out in a manner (block 57) similar to the signal 
update: 

H k = H k + XF *i(g-Htr M ) (s) 

l+XF k :,F k#l 

The same X as before is used in the above equation, which 
is now updated (blocks 58 and 59) through the following 
sequence: impose frequency domain constraints on H k , take 
the inverse Fourier transform of H k , constrain the result 
to be real and positive inside the selected spatial extent 
and zero outside, and normalize it to obtain h k+l , then take 
the Fourier transform of h^, to obtain H^,, where 
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normalization refers to dividing all the elements of h k+1 by 
their sum so that they add up to unity- In some cases, 
such as those encountered in 3-D optical microscopy, 
observed attenuation can be incorporated here by ensuring 
that the impulse response function sums up to an 
appropriate positive number representing an attenuation 
factor which is smaller than unity. The process then 
returns to blocks 53 and 54 to update f k+i . 

A decision is then made at 55 whether to continue 
the iterations. This may simply involve determining 
whether a preselected number "n" of iterations have been 
completed; alternatively, if successive updates of h and f 
do not change by more than prespecified limits, e.g., 



J£*iil£*Ua, and, 



jjW^j ±fL (9) 



where tx x and a 2 are preselected numbers, then the iteration 
continues. When the iterations are completed, f k+l 
constitutes the estimate I and h^, constitutes the estimate 

A variation on this process sequence is shown in 
Fig. 4B. In this alternative sequence, H k + X Is determined 
first at blocks 57, 58 and 59, and then the process 
proceeds to determine f k+l and F k+l at blocks 53, 54 and 56. 
If the selected number of iterations has been completed, as 
determined at block 55, the last updated estimate of f , f a , 
is used as the final estimate i ; if not, the process 
continues through blocks 57, 58 and 59. 

Fig. 5 schematically illustrates the result of 
deconvolution for an exemplary two-dimensional signal which 
represents a two dimensional image. A two-dimensional set 
of signal data g, indicated generally within the dashed 
line 60 of Fig- 5, has a size Nl x N2. This signal data 
set may comprise, for example, a signal corresponding to 
two-dimensional set of pixels which constitute a single 
frame of data from a confocal microscope taken at one focal 
plane, with the value of the signal g at each pixel being 
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the detected light intensity at that pixel. The signal g 
is deconvolved by a point spread function (PSF) h, 
indicated within the line labeled 61 in Fig. 5, which has a 
size Ml x M2, which generally is smaller than the set of 
signal data g. The result of the deconvolution operation, 
to provide the estimate f , is shown within the line 62 of 
Fig. 5, and will be of a size (Nl-Ml+l) x (N2-M2+1) . This 
result will also apply to single dimensional data or to 
three dimensional and higher dimensional data. 

Fig. 6 illustrates the operation of the overlap- 
save method of partitioning the set of signal data for a 
two-dimensional signal. For example, the set of signal 
data 60 of Fig. 5 may be divided into four overlapping 
regions, each one of which is bounded in part by the outer 
boundary 60, the four regions also being bounded by the 
lines 66-69, lines 68-66, lines 68-67, and lines 69-67, 
respectively, if the set of signal data 60 is larger, and 
smaller regions are desired, more than four regions can be 
created. Each of the four regions illustrated in Fig. 6 is 
deconvolved separately in accordance with the process above 
to obtain estimates of the actual signal f within each 
racrion. The resulting regions for which the estimated 
sienalsf are obtained are illustrated by the shaded 
rectangles in Fig. 6. The section of signal data which is 
used to obtain the signal estimate within the shaded 
rectangle is the dotted rectangle which surrounds the 
shaded rectangle. This is similar to the relationship 
between the initial set of signal data 60 and the estimated 
signal data 62 of Fig. 5. With proper choice of the 
regions, the junction between the regions covered by the 
estimated signal is seamless, as indicated in Fig. 6. 
Again, the procedure can be extended to one-dimensional 
data or to three or more dimensional data. 

As indicated above, the parameter X is determined 
based on the noise variance in the signal data g. While X 
will generally be chosen based on the considerations 
discussed above for a particular instrument to which the 
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present invention is applied and the actual noise level, 
the values of X can be empirically determined based on 
experience with a particular instrument. In practice , it 
is typically necessary to scale down the PSF by a factor of 
10 to avoid numerical inaccuracies and then rescale the 
result appropriately to account for this scaling. 
Therefore, the X values used in the discussion below are 
scaled up by 100. Accordingly, the variable sf (scale 
factor) in the computer program listing set forth further 
below is set to 10. For confocal microscopy involving 
three-dimensional signal data which is imaged and 
displayed, the values for X set forth in Table 1 below have 
been found to provide satisfactory performance based on the 
preferred criteria of the user for the visual image 
quality. 



Table 1 



Visual Quality 


X Range 


Low blur and low noise 


5000 - 20000 


Low blur and medium noise 


1000 - 5000 


Low blur and high nnise 


100 - 1000 


Medium blur and low noise 


2000 - 5000 


Medium blur and medium noise 


500 - 2000 


Medium blur and high noise 


100 - 500 


High blur and low noise 


800 - 2000 


High blur and medium noise 


100 - 800 


High blur and high noise 


10 - 100 



As used in Table 1 above, high, medium, and low 
noise levels correspond to signal-to-noise-ratios less than 
15dB, between 15dB and 30dB, and greater than 30dB, 
respectively. As a further example, it has been found that 
most images obtained utilizing the Odyssey* confocal 
microscope system after 3 2 frame averages fall within the 
medium noise category. Thus, a value of X in the range of 
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1000 is typically satisfactory for most applications with 
such an instrument* 

The application of the invention to a simplified 
one-dimensional simulated signal for purposes of 
exemplifying the invention may be illustrated with 
reference to the graphs in Figs. 7A-F. The synthesized 
one-dimensional signal f (the ideal or actual signal) is 
shown in Fig. 7A. It was convolved with the Gaussian- 
shaped impulse response of Fig. 7E (the actual response 
function h) , and zero-mean white noise was then added to 
obtain the degraded signal g shown in Fig. 7B. Although 
the signal g is represented as a function of pixel position 
in Fig. 7B, it could also be a function of any other one- 
dimensional variable, e.g., time. The severity of the 
degradation varies with the spread of the impulse response 
and the noise variance. The Gaussian function of Fig. 7E 
has a variance of 2, and the signal-to-noise ratio in the 
signal g is 3 0dB. The power spectrum of the signal g is 
shown in Fig. 7C. The restored estimated signal £ is shown 
in Fig. 7D, which is seen to be a much closer approximation 
of the input signal f of Fig. 7A than the distorted signal 
g of Fig. 7B. The actual point spieau function (PSF) of 
Fig. 7E may be compared with the estimated point spread 
function shown in Fig. 7F. 

To perform this restoration in accordance with 
the invention as described above, the observed signal data 
g was used as an initial estimate of the signal f and a 
Dirac delta function was used as the original estimate of 
the PSF h. The spatial extent of the PSF was assumed to be 
unknown and was overestimated as 21 pixels wide, even 
though the true PSF was only 13 pixels wide (as shown in 
Fig. 7E) . The spatial extent may be estimated from the 
spread of the transition portion of the peaks shown in the 
signal plotted in Fig. 7B. This spread is seen to be about 
seven or eight pixels, and it is generally preferable to 
overestimate the spatial extent by 200% to 300%. 
Consequently, 21 pixels were chosen. The noise variance 
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was estimated using the portion of the signal g of Fig- 7B 
lying between pixels 10 and 20 (where there was no apparent 
signal data component) • The bandlimit parameter was 
estimated from the power spectrum G of the signal g as 
plotted on a linear scale in Fig. 7C. A bandlimit of 29 
was selected as being the upper frequency limit of useful 
(non-noise) information and was used for this restoration. 
It is seen from Fig. 7C that the magnitude spectrum tends 
to decline from frequency 29 up, with no distinct peaks 
which would indicate the presence of signal information at 
frequencies above 29. Similar estimations may be carried 
out on two- and three-dimensional (or higher dimensional) 
data by plotting the frequency spectrum obtained from 
signal data taken along each orthogonal axis. 

Using the foregoing simulated data, the effects 
of various changes in system conditions on the process can 
be observed. Figs. 8A-8D show the effect of increasing 
noise variance (corresponding to decreases in the SNR from 
40dB to lOdB) in the signal estimate 2. Generally, with 
decreasing SNR, the fidelity of the restorations decreases 
in an expected manner. 

If the PSF variance increases, the fidelity of 
the signal estimate also declines if the spatial extent of 
the PSF is also increased correspondingly. A similar 
decline in accuracy occurs as the PSF is made more 
asymmetric. 

The effect of the bandlimit selection on the 
signal estimate £ is illustrated in the graphs of Figs. 9A- 
9E, which show the results of the estimations $ for 
bandlimits selected at frequencies of 18, 24, 27, 36 and 
45, respectively. The bandlimit used to obtain the 
estimate of Fig. 7D was 29. These simulated results 
indicate that for this signal data, overestimation of the 
bandlimit causes more distortions than underestimating the 
bandlimit in the neighborhood of the correct estimate 
(where 29 is assumed to be the correct estimate) . 
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Generally, it is not as necessary to accurately 
estimate the spatial extent of the PSF as the bandlimit for 
the signal g. In general, an accurate bandlimit estimation 
can automatically estimate the extent of the PSF. Even 
when the PSF extent is deliberately overestimated r a 
correct bandlimit estimate determined from the signal data 
g results in either zeros or insignificantly small values 
beyond the actual extent of the PSF. It should be noted 
that the spatial extent of the PSF must always be 
overestimated in order to take advantage of the influence 
of the bandlimit on the PSF spatial extent. However, the 
true spatial extent can be obtained by initially processing 
a small region or volume of data, and the estimate can then 
be used for the subsequent calculations. 

It is well known in the restoration of 
bandlimited signals that the mean-squared-error decreases 
with iterations up to a certain limit and then starts 
increasing again. By appropriately choosing the bandlimit 
and the noise variance in the present invention as 
discussed above, one can converge in accordance with the 
procedure of the present invention in a fewer number of 
itciation=* aL the expense of a small increase in the noise 
amplification in the final reconstruction. 

In general, for signal-to-noise ratios 
encountered in practice (typically less than 3 0dB) , a PSF 
obtained using geometrical considerations (obtained with 
bandlimit and noise variance parameters determined from an 
analysis of the actual signal data) performs as well as one 
obtained with more rigorous considerations, and that even 
if the exact PSF cannot be reconstructed, if the PSF is 
sufficiently "close 11 , reconstruction of the signal data can 
be obtained which is comparable to that obtained with exact 
prior knowledge of the PSF. 

The simulation described above was carried out 
with one dimensional data. The present invention is also 
applicable to higher spatial dimensions. For example, in 
optical microscopy, the blurring function h is inherently 
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three-dimensional for the reasons described above. For the 
deblurring of signals corresponding to two-dimensional 
images, the following procedures apply: 

(a) Bandlimit netermination. The 2-D frequency 
spectrum may be analyzed along the u (reciprocal of x) and 
v (reciprocal of y) directions, since these directions are 
mutually orthogonal. The bandlimits are determined in 
these orthogonal directions in the same manner as in the 
one-dimensional case (e.g., by display of signal data and 
selection of parameters) and, thus, a rectangular frequency 
bound is established. 

(b) Noise Var iance Determination. Noise 
variance can be determined by choosing for examination a 
rectangular region where the image is visually determined 
to be uniform in intensity. If there is more than one such 
region available, the maximum variance in any of these 
regions is taken as the estimated noise variance. This 
provides a good upper bound for the noise variance. 

There are various situations where it is desired 
to restore a blurry two-dimensional image in optical 
microscopy. Because the point spread function in such 
cases is three-dimensional, it has been suggested in the 
prior art that a projected 2-D point spread function may be 
inadequate. Although this may be true in general, 
substantial improvements in individual 2-D images can be 
obtained using the procedure of the present invention 
followed by adaptive smoothing to suppress noise in uniform 
intensity regions. 

The present invention can also be extended to 
three-dimensional spatial variables, such as in 3-D optical 
microscopy, or higher dimensions. As in the 2-D case, the 
bandlimit is determined for signal data along mutually 
orthogonal directions, e.g., by display of the signal along 
such directions. The procedure is adequate, and may be 
modified if it is determined that a more accurate bandlimit 
is needed. If cones of missing frequencies do not permit a 
cuboidal shape of support, as in widefield microscopy, the 
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procedure may be adapted as required to constrain the 
bandlimits. 

The present invention may be used with a confocal 
microscope system of Fig. 2 for 3-D fluorescence 
microscopy. The procedure is applicable for magnitudes of 
specimens which are obtained either by staining the 
specimen (for transmitted light) or by reflectivity (for 
reflected light) . In transmitted microscopy, the image 
intensities are inverted and the procedure is then applied 
as discussed above. The present invention provides 
accurate estimations of the actual signal data which 
compare well with reconstructions obtained using the more 
laborious empirically determined PSFs for thin specimens. 
Moreover, the present invention provides much improved 
resolution for thick specimens for which, as a practical 
matter, the PSFs cannot be empirically determined. Under 
such conditions, if the prior approaches are used, either 
the data cannot be processed or it is processed with 
inaccurate PSFs and results are obtained which are not much 
improved from the original data. Traditionally, it has not 
been possible to deconvolve thick specimens (as can be done 
with thin specimens), including individual cells, for 
several reasons. First, it was not possible to accurately 
determine empirically the blurring function for a given 
specimen because of several factors which can influence the 
PSF of the observed data, including a refractive index 
mismatch, thickness of the specimen and the attenuation, 
and the wavelength of excitation. Secondly, the blur in a 
thick specimen is spatially variant, which does not match 
the assumptions in most of the known deconvolution 
procedures. Thirdly, the bandlimits on the PSF do change 
because of spatial variance of the PSF, and, therefore, 
theoretically obtained bandlimits are not appropriate for 
thick specimens. 

In the present invention, all of the constraints 
that are applied are derived from the actual signal data g, 
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and 'thus it is possible to apply the present invention to 
thick specimens. 

There have been attempts to describe no-neighbor 
restoration of images which assume that either neighboring 
5 sections do not change in intensity or change linearly in 

intensity , and with the further assumption that the PSF is 
axially symmetric. The blind deconvolution process of the 
present invention can be applied to deblur 2-D images (ray 
summed) to expedite such restoration. 

10 If images have a signal to noise ratio (SNR) less 

than 3 0dB, simulation results have shown that the PSF 
obtained from geometric considerations in accordance with 
the present invention perform comparably to the more 
precisely obtained PSF using other techniques. Thus, for 

15 most practical applications, it is sufficient to use PSFs 

which are close to the true PSFs. If some prior knowledge 
about the PSF is available from empirical studies, the 
empirically determined PSF can be incorporated as the 
starting estimate in the present invention and can be 

20 recursively refined to obtain a more accurate PSF and a 

somewhat improved restoration. Where the empirically 
determined PSF is not available, the Dirac delta function 
is preferred as the starting estimate of the PSF, although 
other functions can be used. For example, a Dirac delta 

25 may be used as the starting estimate of H in the frequency 

domain. 

The buildup of noise in restored images is 
inevitable when the main goal of the restoration is to 
obtain improved resolution. Therefore, the resolution 
30 procedure should optimize both resolution and noise 

enhancement. Generally, it is possible to get improved 
resolution at the risk of slightly higher noise 
amplification. If desired, detail preserving smoothing can 
be preformed. 

35 As noted above, where spatially varying blurs are 

present, the data can be sectioned into smaller segments, 
and for each smaller segment, blur parameters can be 
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individually determined. The procedures of the invention 
can, therefore, be applied when the blur in images varies 
slowly. Such slowly varying degradations are typical of X- 
ray spectroscopy (with increasing keV) , two-dimensional gel 
electrophoresis (with increasing field strength) and in 3-D 
optical microscopy (with increasing axial depth) . 

Generally, the procedures of the present 
invention should not be used where there are abrupt changes 
in the blurring function, unless the response in such 
spatial regions can be determined in other ways, and the 
remainder of the image away from the abrupt change can be 
resolved by segmenting the data as discussed above. 

The present invention thus provides rapid and 
accurate true signal estimation for spatial variables, 
including l-D, 2-D and 3-D signals, allows computation of 
all the necessary parameters from the observed data, and 
typically yields acceptable restorations in very few 
iterations, typically five to ten. The PSF bandlimits can 
be determined in accordance with the invention from the 
power spectrum of the actual signal data. Both the PSF and 
the signal are simultaneously restored. 

An exemplary computer program to be carried out 
within the signal processing system 23, e.g., the 
processing computer 38, to provide the signal processing 
operation of Fig. 4, is set forth below. This processing 
accepts the signal data g which has been stored in digital 
form in the memory 24 and provides an outlet signal, which 
is the estimate £ , to a screen display (such as the video 
screen 4 0) , to provide a two dimensional image of the 
estimated signal. The signal process operates using 
available functional modules, e.g., for performing fast 
Fourier transforms (FFT) and inverse transforms. 
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Copyright © 1994 Noran Instruments, Inc. 
Gopal B. Avinash 

/* 
* 

* This program is based on the blind deconvolution 
algorithm 

of G- Avinash ( 1993 ) . 

* It uses SGI f ft library 
* 

*/ 

# include <stdio.h> 
/ include <math.h> 

/include »noran_data.c w /* Image files reside here 

*/ 

#define MAXARRAY 256 

/define PI 3.1415926 

/define psfsize 13 

/define psfsizez 15 

/define psfsize2 (int)psf size/2 

/define asizeS ( int) psf sizez/2 

/define asize6 psfsizez 

/define freq_xnult 10. 

/define zarray 64 /* Max. number of images 

*/ 

/define SWAP(a,b) tempr«(a) ; (a)=(b) ; (b)=tempr 
float psf [psfsize] [psfsize] [psfsizez] ; 
typedef struct 
{ 

float real,imag; 
} COMPLEX; 
struct { 

int x_origl; 

int y_origl; 

int z_origl; 

int x_orig2; 

int y_orig2; 
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int z_orig2 ; 
} orig_cord; 

COMPLEX *compsf , *compsf ptr , *cominput , *cominputptr ; 
COMPLEX *comoutput , *comoutputptr , * workspace ; 
float * input 1,* input Iptr; 
float sine_array[1024] ; 

int asize2 , asizel , asize, maximum, minimum, k, final ; 
float input_intensity, input_intensityl; 
float lambda=l; 

FILE *fopen(const char *, const char *) , *fp, *f ilein, 
♦fileout; 

float sd,sf,rsf; 

float psf_sum, snr ; 

int 

psf size3 , psf size4 , psf size5 , asize3 , asize4 , x_edge , y_edge , z_ed 
ge; 

int x_size,y_size, z_size,xy_size,ym,yp,xm,xp; 

unsigned char *if ilel, *if ilelptr; 

float gauss_fact[25] ; 

float maxout2 , minout2 , maxout , minout ; 

float x_sum,y_sum, z_sum,y_suml; 

int invert, dat as ize; 

int x_ext,y_ext; 

int imageslices , extra_slices ; 

float alpha; 

main (int argc, char **argv) 

{ 

int x,y,i,j,m,n; 
int fdi,countl,count2; 
if (argc <2) 
{ 

fprintf (stderr, "Usage: sbii-3dd imageslices 
non-invert=0/ invert =l\n w ) ; 
exit(O) ; 

} 

sscanf (argv[l] , "%d M , &imageslices) ; 
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sscanf (argv[2] , w %d", & invert) ; /* Used for images with 
bright 

backgrounds */ 
printf (" Copyright (c) I994\n") ; 
printf (" Gopal B. Avinash\n") ; 

printf (" NORAN INSTRUMENTS # Inc.\n w )? 

printf ( w All Rights Reserved\n w ) ; 

printf ( n \n") ; 

sf - 10-; /* scale factor = 10.0; */ 
rsf - l./sf; 
printf ("\n") ; 

printf ("Starting coordinates (x y z)\n"); 

scanf ( w %d %d %d w , &orig_cord.x_origl, &orig_cord.y_origl, 

fiorig^cord. z_ origl) ; 
printf (»\n M ) ; 
printf ( M x_size=\n w ) ; 
scanf ( M %d w , &x_size) ; 
printf ( M y_size=\n H ) ; 
scanf ( w %d", &y_size) ; 
printf ("z_size=\n M ) ; 
scanf ( n %d w / &z_size) ; 
datasize - x size*y_size*7 - size; 
alpha - l./ (float) datasize; 

inputl - (float *) malloc(datasize*sizeof (float) ) ; 

xy_size =x_size*y_size; 

x_edge = x_size-psf size2 ; 

y_edge » y_size-psf size2 ; 

z_edge = z_size-asize5; 

asize4 = psfsize2; 

asize3 = psfsize; 

orig_cord.x_orig2 - orig_cord.x_origl+ x_size; 
orig_cord.y_orig2 ■» orig_cord.y_origl+ y_size; 
orig_cord.z_orig2 = orig_cord. z_origl+ z_size; 
x_ext * orig_cord.x_orig2 - orig_cord.x_origl; 
y_ext = orig_cord.y_orig2 - orig_cord.y_origl; 
for ( i=or ig_cord . z_or igl ; i<or ig_cord . z_orig2 ; i++) { 
inputlptr = input l+(i-orig_cord.z_or igl )*xy_size; 
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read_image (image_name (i) , i) ; 

} 

printf { "Maximum Value of the input=%f \n" , maxout2) ; 
printf ("Minimum Value of the input=%f \n" ,minout2) ; 
asizel = x — size; 
asize2 * y_size; 

printf ("noise floor value of images^ \n") ; 
scanf ("%f ", &snr) ; 
maxout2 -« snr; 
inputlptr = inputl ; 
for(i=0;i<imageslices;i++) { 
for (y=0;y<y_size;y++) { 

f or (x=0;x<x_size;x++, inputlptr++) { 
*inputlptr = (* input Iptr-snr ) ; 
if (*inputlptr <0.) *inputlptr = o.; 
if ( maxout2< *inputlptr) maxout2= (*inputlptr) ; 
if( minout2> *inputlptr)minout2=(*inputlptr) ; 
} 

) 

> 

printf ( "Maximum Value of the corrected 
input=%f \n",maxout2) ; 

printf ("Minimum Value of the corrected 
input=%f \n",minout2) ; 

/it******************************************** j 

printf ( "Enter the number of iterations wanted\n") ; 

scanf ("%d",&final) ; 

for ( i=0 , x=l ; i<=asize4 ; i++ , x++) 

gauss_fact[i] = exp(- (float) (x*x/asize3) ) ; 
inputlptr = inputl; 
for(i=0;i<z_size;i++) { 
for (y=0;y<y_size;y++) { 

for (x=0;x<x_size;x++, inputlptr++) { 
if ( (x>=x_edge) j | (y>=y_edge) } j (y<asize4) \ J (x<asize4) ) 

{ 

if ( (x>=x_edge) && (y>=asize4) && (y<y_edge) ) 
♦inputlptr » * (inputlptr-x+x_edge-l)+ 
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(* (inputlptr-x+asize4 ) - * ( inputlptr-x+x_edge-l) ) 
* (x-x_edge+l) / (float) asize3 ; 
else if ((x<asize4)&&(y>=asize4)&&(y<y_edge)) 

♦inputlptr = *( input lptr-x+x_edge-l)+ (-* (input lptr-x 
+x_edge-l) + * (inputlptr-x+asize4) ) * (asize4+x) / 
(float) asize3 ; 
else if ( (y>=y_edge) && (x>=asize4) && (x<x_edge) ) 
♦inputlptr = * ( inputlptr- (y-y_edge+l) *x_size) + 
(- * ( inputlptr- (y-y_edge+l ) *x_size) + 
*( inputlptr- (y-asize4) *x_size) ) * (y-y_edge+l) / (float) asize3 
else if ( (y<asize4 ) && (x>=asize4 ) && (x<x_edge) ) 
♦inputlptr = * ( inputlptr- (y-y_edge+l) *x_size) + 
( *( inputlptr- (y-asize4 ) *x_size) - 
*( inputlptr- (y-y_edge+l)*x_size) ) *(asize4+y) / (float) asize3 
else if ( (x<asize4)&&(y<asize4) ) *inputlptr = 
* (inputlptr- (x-asize4) - (y-asize4) *x_size) * 
gauss_fact[asize4-x]*gauss_fact[asize4-y] ; 
else if ((x>=x_edge)&&(y>=y_edge) ) *inputlptr - 
* (inputlptr- (x-x_edge+l) - (y-y_edge+l) *x_size) * 
gauss_f act[asize4-x_size+x] *gauss_f act (asize4-y_size+y ] ; 
else if ( (x>=x_edge)&&(y<asize4) ) *inputlptr = 
* ( inputlptr- (x-x_edge+l ) - (y-asize4 ) *x_size) * 

gauss_fact[asize4-x_size+x]*gauss_fact[asize4-y]; 

else if ((x<asize4)&&(y>=y_edge) ) *inputlptr = 

*( inputlptr- (x-asize4 ) - (y-y_edge+l) *x_size) *gauss_ 
fact(asize4-x]*gauss_fact[asize4-y_size+yj ; 



> 



} 

} 

} 

inputlptr=inputl ; 

extra_slices - z_size-imageslices; 
for ( i=0 ; i<z_size ; i++) { 
for (y=0 ;y<y_size ;y++) { 

for (x=0;x<x_size;x++, inputlptr++) { 

if (i>=z_edge) { 

♦inputlptr - *( inputlptr- (i-z_edge+l)*xy_size)+ 



WO 96/10794 



PCT/US95/12259 



- 35 - 

(* (inputlptr-(i-asizeS) *xy_size) - 

* (inputlptr-(i-z_ edge+1) *xy_size) ) * (i-z_edge+l) / 
( (float) asize6) ; 

} 

5 else if (i<asize5) { 

*inputlptr = * ( inputlptr- ( i-z_edge+l) *xy_size) + 
( *{ inputlptr- (i-asize5) *xy_size)- 

* (inputlptr- (i-z_ edge+1) *xy_size) ) * (i+asize5) / 
( (float) asize6) ; 

10 } 

} 

> 

} 

comoutput = (COMPLEX *) malloc (datasize*sizeof (COMPLEX) ) ; 
1 5 input lptr=inputl ; 

comoutputptr=comoutput ; 
f or (i=0;i<z_size;i++) { 
for (y=0 ;y<y_size;y++) { 

for (x=0;x<x_size;x++ , inputlptr++, comoutputptr ++ ) { 
20 comoutputptr->real = *inputlptr; 

> 

} 

} 

/* Low-pass filter the edges of image with 3X3 mask */ 
2 5 input lpt r« input 1 ; 

comoutputptr=comoutput ; 
for ( i=0 ; i<z_size ; i++) { 
f or (y=0;y<y_size;y++) { 

for (x=0;x<x_size;x++, inputlptr++, comoutputptr++) { 
30 if ( (x>=x_size-asize4-l) | { (y>=y_size-asize4-l) } | 

(y<=asize4) { | (x<=asize4)) 

{ 

if (x==0) xm=x_size-l ;else xn « x-1; 
if (y— 0)ym-y_size-l;else ym = y-1; 
35 if (x»»x_size-l) xp«0 ; else xp = x+1; 

if (y=y_size-l) yp=0 ; else yp = y+1; 
comoutputptr->real = ( * ( inputlptr- ( x-xm) ) + 
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* (inputlptr (x-xm) -(y-yp) *x_size)+ 

* ( inputlptr- (x-xm) - (y-ym) *x_size) + 
♦inputlptr + * (inputlptr- (y-yp) *x_size) + 
* ( inputlptr- (y-ym) *x_ size) + 

5 * ( input Ip tr- ( x-xp ) ) + 

* ( inputlptr- (x-xp) - (y-yp) *x_size) + 
♦(inputlptr- (x-xp) -(y-ym) *x_size) ) /9. ; 

} 

else comoutputptr->real « *inputlptr; 
10 } 
} 

} 

free ( input 1) ; 

15 cominput = (COMPLEX *) malloc (x_size*y_size*z_size 

*sizeof (COMPLEX) ) ; 
compsf = (COMPLEX *) malloc (x_size*y_size*z_size 

*sizeof (COMPLEX) ) ; 
workspace « (COMPLEX *) 
20 malloc ( (x_size+y_size+z_size+45) *2 

*sizeof (COMPLEX) ) ; 
cf ft3di (x_size,y_size, z_sizc, workspace) ; 
compsf ptr = compsf; 
cominputptr = cominput; 
25 comoutputptr = comoutput; 

for (i-0;i<z_size;i++) { 
for (y=0;y<y_size;++y) { 

for ( x=0 ; x<x_s i z e ; ++x , ++comps f ptr , ++cominputptr , 
++comou tputptr ) 

30 { 

cominputptr->real « comoutputptr->real; 
compsfptr->real» rsf ; /* Take 3-D fft of psf */ 
compsf ptr->imag=cominputptr->imag=comoutputptr->imag=0; 
} 

35 } 
} 

maxout2=0. ; 
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minout2 =» 10000.; 
comoutputptr*comoutput ; 
for(i=0;i<z_size;i++) { 
for(y=0;y<y_size;y++) { 

for (x-0;x<x_size;x++, comoutputptr++) { 
if (comoutputptr->real<0) comoutputptr->real =0; 
if ( comoutputptr->real>maxout2 ) maxout2= comoutputptr- 
real; 

if ( comoutputptr->real<minout2 ) minout2* comoutputptr- 
real; 

> 

} 

> 

pr intf ( "Maximum Value of the input=%f \n" ,maxout2) ; 
printf ("Minimum Value of the input=%f \n" ,minout2) ; 
for ( i=0 ; i<z_s ize ; i++ ) 

wr ite_image ( input_name ( i ) , i , 1 ) ; 
/* -1 is for the forward fft */ 

cfft3d(-l,x_size,y_size,z_size,cominput,x_size,y_size, 
workspace) ; 

Printf ("The band-limit support parameter is very 
important 
\n"); 

printf ("for the convergence of the algorithm. Choose the 
\n") ; 

printf ("value so that you overestimate it (in your 
opinion) 
\n") ; 
printf ("\n") ; 
cominputptr = cominput; 
comoutputptr = comoutput; 
x_sum =0; 

for ( i»0 ; i<z_size ; i++) { 
for (y=0;y<y_size;++y) { 

for ( x=0 ; x<x_size ; ++x, ++cominputptr , ++comoutputptr ) 

{ 

comoutputptr->real = cominputptr->real; 



WO 96/10794 



PCT/US95/12259 



- 38 - 

comoutputptr->imag - cominputptr->imag ; 
if ( <x<x_size/2 ) && (i=0) && (y=»0) ) 
printf ("%d 

%f \n" # x r hypot(comoutputptr->real f comoutputptr->imag) ) ; 
} 

} 

> 

printf ("band-limit support for x= \n w ) ; 
scanf («%d",&psfsize3) ; 
cominputptr = cominput; 
comoutputptr « comoutput; 
y_sum =0; 

for (i=0;i<z_size;i++) { 
for (y=0;y<y_size;++y) { 

for ( x=0 ; x<x_size; ++x , ++cominputptr , ++comoutputptr ) 

{ 

if ( (y<y_size/2) && (x=0) && (i==0) ) 
printf ("%d 

%f \n" # y # hypot(comoutputptr->real # comoutputptr->imag) ) ; 
} 

} 

* 

printf ("band-limit support for y= \n M ) ; 
scanf («%d",&psfsize4) ; 
cominputptr = cominput; 
comoutputptr = comoutput; 
z_sum =0; 

for ( i=0 ; i<z_size ; i++) { 
for (y=0;y<y_size;++y) { 

for (x=0 ; x<x_s ize ; ++x , ++cominputptr , ++comoutputptr ) 

{ 

if ( (i<z_size/2 ) && (x=0) && (y«=0) ) 
printf ("%d 

%f \n" , i # hypot (comoutputptr->real # comoutputptr->imag) ) ; 
} 
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printf ("band-limit support for z= \n n ) ; 
scanf ( w %d M , &psf sizes ) ; 
/********************************************************* / 

deconv ( ) ; 

} 

/******************* En<a of Main ********************/ 

read_image(char *buf , int i) 
{ 

unsigned char *if ilel, low_byte,hi_byte; 
unsigned char *ifilelptr; 
int col_len,row_len; 
int fdi,x r y; 

fdi = open(buf,0); 
if (fdi<0){ 

pr int f ("error in opening the f ile\n M ) ; 
exit(l); 

> 

asizel=x_size; 
asize2=y_size ; 

if ilel» (unsigned char *)malloc(asizel*asize2*sizeof 

(unsigned char) ) ; 
read(fdi, if ilel,asizel*asize2*sizeof (unsigned char)) ; 
if ilelptr=if ilel ; 
if ( inverts 1) { 

f or (y=0;y<asize2;y++) { 

for (x=0;x<asizel;x++, ifilelptr++) { 
if ( (y>-orig_cord-y_origl) && 
(y< orig_cord.y_orig2) && 
(x>=orig_cord.x_origl) && 
( x< or ig_cord . x_or ig2 ) ) { 
*inputlptr - 255- (float) *ifilelptr; 
if( maxout2< *inputlptr)maxout2« *inputlptr ; 
if( minout2> * input lptr)minout2= *inputlptr ; 
inputlptr++; 
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} 

} 

} 

} 

else{ 

for (y=0;y<asize2 ;y++) { 

for(x=0;x<asizel;x++,ifilelptr++) { 

if ( (y>«orig_cord.y_origl) &S 
(y< orig_cord . y_orig2 ) && 
(x>=orig_cord.x_origl) S& 
( x< or ig_cord . x_or ig2 ) ) { 
♦inputlptr = (float) *ifilelptr; 

if( maxout2< * input lptr ) maxout2= * input lptr ; 
if( minout2> *inputlptr)minout2= *inputlptr ; 
inputlptr++; 

} 

} 

} 

> 

free(ifilel) ; 
close (fdi) ; 



write_image(char *bufl, int i, int spec) 

z*******************************************************/ 

{ 

register int fd,x,y,asize; 
unsigned char *buf , *buf ptr ; 
float *ipixptr,pix_int; 

asize = (asizel>asize2) ? asizel: asize2; 
buf= (unsigned char *)malloc(asize*asize*sizeof 

(unsigned char) ) ; 
f d=cr eat (buf 1,0755) ; 
if (fd<0) { 

printf ("error in creating the file\n") ; 
exit(i) ; 
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buf ptr^buf ; 

comoutputptr=comoutput+i*xy_size ; 
if (spec«=l) { 

for (y«0;y<asize;y++) { 

for (x=0;x<asize;x++,bufptr++, comoutputptr++) { 
if(maxout2 = 0) { 
*bufptr =0; 
continue; 

} 

/* to truncate the edges */ 
else 

if ((x>=x_size-(psfsize2)) j { (y>=y_size- (psf size2) ) | | 
(x<=psfsize2) ! j (y<=psfsize2) { } 
(i>=z_size-(asize5) ) | \ (i<=asize5) ) 
♦bufptr =0; 
else if ( (x<x_ext) &&(y<y_ext) ) { 
if ((pix_int= 
(comoutputptr->real*255./maxout2) ) >255. ) { 
pix_int = 255*; 

} 

*bufptr= (unsigned char) pix_int ; 

ii-(invert=l) *bufptr= (unsigned char) 255- 

*bufptr; 

} 

else { 

*bufptr=0; 

} 

} 

> 

} 

else { 

for (y»0;y<y_size;y++) { 

for (x=0 ;x<x_size;x++ , buf ptr++, comoutputptr++) { 
if (maxout2 =» 0) { 
*bufptr =0; 
continue; 

} 
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/* to truncate the edges */ 
else 

if ( (x>=x_size-(psfsize2) ) | \ (y>=y_size- (psf size2) ) | | 
(x<=psfsize2) ! ! (y<~psfsize2) | | 
5 (i>=z_size-(asize5) ) | | (i<=asize5) ) 

*bufptr =» 0; 
else if ( (x<x_ext ) && (y<y_ext) ) { 
*buf ptr= (unsigned 
char) ((float) (comoutputptr->real-minout2) *255. / 
10 (maxout2-minout2) ) ; 

if (invert=l) *bufptr = 255- *bufptr; 

> 

else { 

*buf ptr=0 ; 

15 } 

} 

> 

} 

write(fd,buf ,asize*asize*sizeof (unsigned char) ) ; 
20 free(buf); 

close (fd) ; 

} 

/★a*******************************************************/ 

write_psf (char *bufl, int i) 
25 ^ **************** **************************************** / 

{ 

register int fd,x,y; 
unsigned char *buf , *bufptr ; 

buf« (unsigned char *)malloc(psf size*psf size*sizeof 
30 (unsigned char) ) ; 

fd=creat(buf 1,0755) ; 
if (fd<0){ 

printf ("error in creating the file\n M ); 
exit(l); 

35 } 

buf ptr=buf ; 

for (y=0;y<psfsize;y++) { 
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for (x==0;x<psfsize;x++,bufptr++) { 
if (maxout2 — 0) { 
♦bufptr * 0; 
continue ; 
} 

*buf ptr» ( unsigned 
char) ( (float) (psf [x] [y] [i]-minout2) *255./ 

(maxout2-minout2) ) ; 
if (invert=»l) *bufptr « 255- *bufptr; 
> 

> 

vrite(fd,buf ,psfsize*psfsize*sizeof (unsigned char)) ; 
free(buf); 

close (fd) ; 

} 

deconv(void) 
{ 

int x,y, i,j , count, countl; 
float lambdal; 
float psf_energy; 
float mag2psf ; 
if (k«=0){ 

lambda= 1000. ; 

printf ("lambda « %f \n»Manbda) ; 
printf ("required lambda - \n") ; 
scanf ("%f",&lambda) ; 

> 

compsfptr = compsf ; 

cominputptr = cominput; 

comoutputptr = comoutput; 

for (i=0;i<z_size;i++) { 
for (y~0;y<y_size;++y) { 
for ( x=0 ; x<x_size ; x++ , compsf ptr++ , cominputptr ++ , comoutputptr 
++) 

{ 

mag2psf » compsfptr->real*compsfptr->real+ 
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comps f p t r - > imag * c omp s f p t r - > imag ; 
comoutputptr->real = 
(comoutputptr->real+lambda* ( (compsfptr->real * 

cominputptr->real) + 

(compsf ptr->imag * cominputptr->imag) ) ) / 

(1+ lambda* mag2psf ) ; 
comoutputptr->imag = 
(comoutputptr->imag+lambda*( (compsf ptr->real * 

cominputptr->imag) - 

(compsfptr->imag * cominputptr->real) ) ) / 

(1+ lambda * mag2psf ) ; 

} 

} 

} 

cf f t3d ( 1 , x_size , y_size , z^size, comoutput , x^size, y_size , works 
pace) ; 

comoutputptr = comoutput; 
for (i=0;i<z_size;i++) { 
for (y=0;y<y_size;y++) { 

for (x=0;x<x_size;x++,comoutputptr++) { 
comoutputptr- >real /= (float) datasize; 
como-Jtputptr->imag /= (float) datasize; 
} 

} 

} 

comoutputptr=comoutput ; 
for ( i»0 ; i<z_size ; i++) { 
for (y=0;y<y_size;y++) { 

for ( x=0 ; x<x_size ; x++ , comoutputptr++) { 
if ((x>=x_edge) { | (y>=y_edge) | | (y<asize4) ■ | (x<asize4)) 

{ 

if ( (x>=x_edge) && (y>=asize4) && (y<y_edge) ) 
comoutputptr->real « 
(comoutpiitptr- (x-x_edge+l) ) ->real+ 

(- (comoutputptr- (x-x_edge+l) ) ->real+ 
(comoutputptr- (x-asize4 ) )->real)* 
(x-x_edge+l) / (float) asize3 ; 
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else if ((x<asize4)S&(y>=asize4)&S(y<y_edge)) 
comoutputptr->real = 
(comoutputptr- (x-x_edge+l) ) ->real+ 

(- (comoutputptr- (x-x_edge+l) ) ->real+ 
(comoutputptr- (x-asize4 ) ) ->real) * 
(asize4+x) / (float) asize3 ; 
else if ( (y>=y_edge) && (x>-asize4) && (x<x_edge) ) 
comoutputptr->real «■ 
(comoutputptr- (y-y_edge+l) *x_size) ->real+ 

(- (comoutputptr- (y-y_edge+l) *x_size) ->real+ 
(comoutputptr- (y-asize4 ) *x_size) ->real) * 
(y-y_edge+i) / (float) asize3 ; 

else if ((y<asize4)&&(x>=asize4)&&(x<x_edge)) 
comoutputptr->real = 
(comoutputptr- (y-y_edge+l) *x_size) ->real + 

( (comoutputptr- (y-asize4 ) *x_size) ->real - 
(comoutputptr- (y-y_edge+l) *x_size) ->real) * 
(asize4+y) / (float) asize3 ; 
else if ((x<asize4)&&(y<asize4))comoutputptr->real 

(comoutputptr- (x-asize4)-(y-asize4)*x_size)->real * 
gauss_ f dct[asi^e4-x] * 

gauss_fact [asize4-y] ; 

else if ( (x>=x_edge) && (y>=y_edge) ) comoutput->real 

(comoutputptr- (x-x_edge+l) - (y-y_edge+l) ) ->real * 
gauss_f act (asize4-x_size+x J * 
gauss_f act[asize4-y_size+y] ; 
else if ((x>=x_edge)&&(y<asize4) ) 
comoutputptr->real= 

(comoutputptr- (x-x_edge+l) - (y-asize4) ) ->real* 
gauss_f act [asize4-x_size+x] * 
gauss_fact[asize4-y] ; 
else if 

( (x<asize4) && (y>=y_size-aslze4) ) comoutputptr->real - 

( comoutputptr- (x-asize4) - (y-y_edge+l) ) ->real * 
gauss_f act[asize4-x] * 
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gauss_f act[asize4-y_size+y] ; 

} 

> 

} 

} 

comoutputptr=comoutput ; 
for ( i-0 ; i<z_size ; i++) { 
for(y«0;y<y_size;y++) { 

for (x-0 ;x<x_size;x++, comoutput:pt:r++) { 
if (i>«z_edge){ 
comout:put:ptr->real « 
(comoutputptr- ( i-z_edge+l) *xy_size) ->real+ 

( (comoutputptr-(i-asize5) *xy_size) ->real - 
(comoutputptr-(i-z_edge+l) *xy_size) ->rcal) * 
(i-z_edge+l) / ( (float) asize6) ; 

} 

else if (i<asize5) { 

comoirtputptr->real » 
(comoutputptr- ( i-z_edge+l) *xy_size) ->real + 

( (comoutputptr-(i-asizeS) *xy_size) ->real - 
(comoutputptr- (i-z_edge+l) *xy_size) ->real) * 
(i+asize5) / ( (float) asiscS) ; 

} 

} 

} 

} 

comoutputptr=comoutput ; 
for(i=0;i<z_size;i++) { 
for (y=o;y<y_size;y++) { 

f or (x=0 ;x<x_size;x++, comoutputpt:r++) { 
if ( (x>=x_size-asize4-l) | J (y>«y_size-asize4-l) ! j (y<=asize4) 
| (x<=asize4) ) 
{ 

if (x=o)xni=x_size-l;else xm = x-1; 
if (y™o)ym-y_size-l;else ym = y-l; 
if (x=x_size-l) xp=0 ; else xp = x+1; 
if (y==sy_size-l)yp=0; else yp = y+1; 
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comout:putptr->iinag = ( ( comoutputptr- (x-xn) ) ->real + 

(comoutputptr- (x-xa) - (y-yp) *x_size) ->real h 
(comoutputptr- (x-xm) - (y-ym) *x_size) ->real+ 
comoutputptr->real+ ( comoutputptr- (y-yp) 

*x_size)->real+ 
(comoutputptr- (y-ym) *x_size) ->real + 
(comoutputptr- (x-xp) ) ->real+ (comoutputp'tr- 
(x-xp) - (y-yp) *x_size) ->real+ 
(comoutputptr- (x-xp) - (y-ym) *x_size) ->real) /9 . ; 



} 

} 

} 

comoutputptr = comoutput ; 
15 for (i=0;i<z_size;i++) { 

for (y=0;y<y_size;y++) { 

for (x=0;x<x_size;x++, comoutputptr++) 

{ 

if ( (x>=x_size-asize4-l) J | (y>=y_size-asize4-l) { | 
20 (y<«asize4) \ \ (x<=asize4) ) 

{ 

comov/tputptr->real = comoutputptr->imag; 

if (comoutputptr->real < 0) comoutputptr->real — 0; 

comoutputptr->imag = 0; 

25 } 

} 

} 

} 

comoutputptr^comoirtpu-t ; 
30 for(i=0;i<z_size;i++) { 

for (y=0;y<y_size;++y) { 

f or (x=0 ;x<x_size;x++,comoutputptr++) { 
if (comoutputpt:r->real<0)comoutputptr->real » 0; 
comoutputptr->imag - 0; 
35 if (maxout2<comoutputptr->real) maxout2 = 

comoutputptr- >real ; 
} 
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} 

} 

printf ( w lteration«%d , maxout»*f\n n , k , maxout2 ) ; 
for ( i=0 ; i<z_size ; i++) 

vr it e_image ( output_name ( i ) , i , 2 ) ; 
prixvtf ("max « %f, min - %f \n» ,maxout2 ,minout2) ; 
if (k«-final)exit(0) ; 
/★A******************************************************/ 
cf f 1 3d ( -1 , x_size , y_s ize , z_size , comoutput , x_size , y_size , 
workspace) ; 

if (k>=10){ 

goto skip_psf ; 

} 

compsfptr = compsf; 
cominputptr = cominput; 
comoutputptr «= comoutput; 
/* Compute the signal */ 
psf_energy = 0; 
for(i=0;i<z_size;i++) { 
for (y=0;y<y_size;++y) { 
for (x=0 ; x<x_si*e ; x-r+ , comoutputptr++ , compsf ptr++ , cominputptr 

++) 

{ 

mag2psf = comoutputptr->real*comoutputptr->real+ 

comoutputptr-> imag * comoutputptr- > imag ; 
if (((X<=psfsize3)&&{y<=psfsize4)&&(i<=psfsize5)) | j 
( (x>=x_size-psf size3) && (y<«psf size4) && 

(i<=psf size5 ) ) j { 
( (x<=psf size3) &6 (y>«y_size-psf size4) && 

(i<=psf sizes)) J S 
( (x>=x_size-psf size3 ) && (y>=y_size-psf size4 ) && 

(i<=psfsize5)) | S 
( (x<=psf size3 ) && (y<=psf size4 ) && ( i>=z_ 

size-psf sizes) ) j j 
( (x>«x_size-psf size3) && (y<=psf size4) && (i>=z_ 
size-psf sizes) ) J J 
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( (x<=psf size3) && (y>=y_size-psfsize4) && (i>=z_ 

size-psfsize5) ) j j 
( (x>=x_size-psf size3 ) && (y>=y_size-psf size4 ) && 

(i>=z_size-psf sizes) ) ) 

{ 

compsfptr->real = (compsf ptr->real+lambda 
* ( ( comoutputptr-> 

real * cominputptr->real) + comou*tputptr->imag * 
cominputptr->imag) ) ) / 

lambda * (mag2psf ) ) ; 
compsf ptr->imag = ( compsf ptr->imag+lambda * ( ( comoutputptr 
->real * cominputptr->imag) - (comoutputptr->imag * 
cominputptr->real) ) ) / 

(1.+ lambda * (mag2psf ) ) ; 

} 

else { 

psf_energy +* 

compsfptr->real*compsfptr->real+compsfptr->imag*compsf p*tr-> 
imag; ; 

compsf ptr->real =* compsf ptr->imag = 0;} 

} 

> 

} 

printf( M psf energy outside «%f \n" ,psf_energy) ; 
/*★**★★**** ******************************* ******* ********* / 

cfft3d(l,x_size,y_size,z_size, compsf ,x_size,y_size,vorkspac 

e); 

compsfptr = compsf; 
for ( i«0 ; i<z_size ; i++) { 
for (y«0;y<y_size;y++) { 

for (x=0;x<x_size;x++, compsf ptr++) { 
compsf ptr->real /= (f loat)datasize; 
compsfptr->imag /=» (float) datasize; 
} 

} 

} 

psf_sum » 0; 
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compsf ptr»compsf ; 
for ( i=0 ; i<z_size ; i++) { 
for(y-0;y<y_size;y++) { 

f or (x«0 ;x<x_size ;x++, compsf ptr++) 

{ 

if ( { (x<=psfsize2)&&(y<»psfsize2)&&(i<=psfsizez) ) | | 
( (x>=x_size-psfsize2) 6& (y<»psf size2) && (i< 

spsfsizez) ) ! ! 
( ( X <«psf size2) && (y>=y_size-psf size2) && (i< 

=psf sizez) ) J ! 
( (x>=x_size-psf size2 ) && (y>=y_size-psf size2 ) && ( i< 

spsfsizez) ) ! ! 
( (x<=psf size2) && (y<=psf size2) *&(i>=z_size 

-psf sizez) ) j | 
( ( x>=x_size-psf size2 ) && (y<=psf size2 ) && ( i>=z_size 

-psf sizez) ) ! ! 
( (x<=psf size2) && (y>=y_size-psf size2) && (i>=z_size 

-psf sizez) ) ! j 
( (x>=x_size-psf size2) && (y>=y_size-psf size2) &&(i>=z 

— size-psf sizez) ) ) 

{ 

if (compsfptr->real<0)compsfptr->real = 0; 
compsfptr->imag =0; 

> 

else{ 

compsfptr->imag = 0;compsf ptr->real = 0;} 

} 

} 

} 

compsf ptr«compsf ; 
for (i=0;i<z_size;i++) { 
for(y=0;y<y_size;y++) { 

for (x=0;x<x_size;x++, compsf ptr++) 

{ 

psf_sum +- compsf ptr->real; 

} 

> 
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} 

printf ("psfsu* « %f \n'\psf_sum*sf ) ; 
compsf ptr=compsf ; 
for(i=0;i<z_size;i++) { 
for (y=0;y<y_size;y++) { 

for ( x=0 ; x<x_size ; x++ , compsf ptr++ ) 

{ 

compsfpt:r->real /=* (psf_sum*sf ) ; 
compsf ptr->imag = 0; 

} 

} 

} 

/* As defined earlier, for the following loop, asizeS 

= psfsizez/2 */ 
compsf ptr = compsf; 
for(i«=0;i<z_size;i++) { 
for (y=0;y<y_size;++y) { 

for (x=0 ;x<x_size ; ++x , ++compsf ptr ) 

{ 

if (((x<=psfsize2)&&(y<=psfsize2)&&(i<=asize5)) \ \ 
((X>=x_size-psfsize2)&&(y<=psfsize2)&&(i<=asize5)) j | 
i Cx<^psfsize2) S&(y>=y_size-psfsize2)&&(i<=asize5) ) | J 
( Cx>=x_size-psfsize2) && (y>«y_size-psf size2) && 
(i<=asize5) ) \ \ 

( (x<=psfsize2) &&(y<=psfsize2) &&(i>=z_size-asize5) ) | | 
( (x>=oe_size-psf size2) &&(y<=psf size2) && (i>=z 
_size-asize5) ) } J 

( (x<=psfsize2) &&(y>«y_size-psfsize2) && (i>=z_size 
-asizeS) ) \ \ 

( (x>=x_size-psf size2) && (y>«y_size-psf size2) && 
(i>*=z_size-asize5) ) ) 
psf [ (x+psfsize2)%x_size] [ (y+psf size2) %y_size] [ (i+asize5) %z_ 
size]= compsf ptr->real; 
} 

} 

} 

f or (i=0 ; i<psf sizez ; i++) 
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write_psf (psf_name (i) , i) ; 
cf f t3d ( -1 , x_size , y_size , z_size , compsf , x_size , y_size , workspa 
ce) ; 

5 skip_psf: 

if (k<«final) { 
k+«l; 
deconv() ; 

} 

10 } 



It is understood that the invention is not 
confined to the particular embodiments set forth herein as 
illustrative, but embraces all such forms thereof as come 
within the scope of the following claims. 
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CLAIMS 



What is claimed is: 

1. A method of reconstructing the output signal 
g from an acquisition system instrument to estimate the 
ideal input signal f to the acquisition system instrument, 
comprising the steps of: 

(a) using the instrument to acquire a 
signal g which is dependent on a variable of one or more 
dimensions, and storing the spatial domain signal g in 
digital form in a memory of a signal processor; 

(b) determining using the signal processor 
the Fourier transform G of the output signal g, and 
selecting frequency bandlimits in each dimension based on 
the range of signal frequency content exhibited by G in 
each dimension; 

(c) selecting an initial response function 
h c for the acquisition system, and determining the Fourier 
transform of h e and applying in the signal processor the 
selected bandlimits to the Fourier transform to provide an 
initial frequency domain estimate h o ; 

(d) then iteratively determining in the 
signal processor estimates of the ideal input signal f and 
of the actual acquisition system instrument response 
function h, where the number of the last iteration is n, 
comprising at the kth iteration the steps of: 

(1) determining an estimator F k in the 

frequency domain as 

F k * Pfc «. *£jg^w 

1 ♦ XH k *H k 

for Jc m o, l , n 

F. - G, 

H, is determined as set forth above, 
H k is the complex conjugate of 1^; 
X is a selected parameter. 
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(2) determining an estimate f k+1 for 
the ideal input signal as the inverse Fourier transform of F k 



3 for all values thereof which are real and positive within a 

4 selected spatial extent and as zero otherwise; 

5 (3) determining the Fourier transform 

6 F k +, of f k>l ; 

7 (4) determining an estimator H k as 

H> - H> ^.--<°-".r.-.> 

1 + XFktiF^j 

8 where F # k+1 is the complex conjugate of F k+1 ; 

9 (e) applying the selected frequency 

10 bandlimits to H k ; and 

11 (f) determining an estimate h k+l of the 

12 acquisition system instrument response function as the 

13 inverse Fourier transform of the bandlimited H k for all 

14 values thereof which are real and positive within a 

15 selected spatial extent and as zero otherwise and 

16 normalized to a selected value; and repeating the steps 

17 until a selected criterion is met and providing the final 

18 estimate f a as an output signal which is an improved 

19 estimate of the ideal input signal f. 

1 2. The method of Claim 1 wherein the initial 

2 estimate ho of the acquisition system instrument response 

3 function is the Dirac delta function. 

1 3. The method of Claim 1 wherein the parameter 

2 X is selected based on the noise variance of the signal g, 

3 with lower values of X selected where high noise variance 

4 is present and higher values of X selected where lower 

5 noise variance is present. 

1 4 # The method of Claim 1 wherein the parameter 

2 X is selected to have a value in the range of 100 to 

3 20,000. 
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5. The method of Claim 1 wherein the parameter 
X is selected to have a value of approximately 1,000. 

6. The method of Claim 1 wherein the output 
signal g of the acquisition system instrument is a function 
of a variable which is at least two dimensional, and 
further including the steps of storing data corresponding 
to the signal data g in a memory of the signal processor, 
subdividing the signal data g into a plurality of sets of 
signal data, the signal data composing each set partially 
overlapping the signal data for each adjacent set, and 
carrying out the method of Claim 1 on each signal data set 
to determine a response function fi and estimated signal 
function i for each signal data set. 

7. The method of Claim 6 wherein the variable 
of which the output signal g is a function is three 
dimensional. 

8. The method of Claim 1 wherein the 
acquisition system instrument is a confocal microscope 
system and wherein the output signal g comprises digital 
data corresponding to light intensity values in a three 
dimensional array corresponding to light received from a 
three dimensional sample examined by the confocal 
microscope system. 

9. The method of Claim 8 including the further 
step of displaying on a video screen to an operator a 
visual image of the estimated ideal input function £ for 
data points corresponding to a single layer in a three 
dimensional array of signal data obtained from the sample 
by the confocal microscope system. 

10. The method of Claim 1 wherein in the step of 
normalizing the response function h^, the response 
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1 function is normalized so as to have a sum of values equal 

2 to one. 

1 11. The method of Claim 1 wherein the criterion 

2 for the last iteration is that a selected number of 

3 iterations is completed which is less than or equal to 10. 

1 12. A method of reconstructing the output signal 

2 g obtained from an acquisition system instrument to 

3 estimate the ideal input signal f to the acquisition system 

4 instrument, comprising the steps of: 

5 (a) using the instrument to acquire a 

6 signal g which is dependent on a variable of one or more 

7 dimensions; 

8 (b) providing the signal g from the 

9 instrument to a signal processor and storing the spatial 

10 domain signal g in digital form in the memory of the signal 

11 processor; 

12 (c) estimating for the signal data g a 

13 selected bandlimit of the significant frequencies in the 

14 signal data g and providing that estimate to the signal 

15 processor; 

16 (d) beginning in the signal processor an 

17 iterative procedure by estimating in the frequency domain 

18 the ideal input signal f using a selected initial estimate 

19 of the impulse response function h of the instrument 

20 transformed to the frequency and the signal data g 

21 transformed to the freguency domain, transforming the 

22 frequency domain estimate to the spatial domain and 

23 applying constraints in the spatial domain to the estimate 

2 4 of f to constrain the estimate of f to be real and positive 

25 within a selected spatial extent and otherwise zero, using 

26 the transform to the frequency domain of the constrained 

27 estimate in the spatial domain of f , the transform to the 

28 frequency domain of the signal g and the previous estimate 

29 of h in the frequency domain to obtain an estimate of h in 

30 the frequency domain which is constrained in the frequency 



WO 96/10794 



PCT7US95/12259 



- 57 - 

domain to be within the selected bandlimit, transforming 
the constrained frequency domain estimate to the spatial 
domain and constraining the estimate to be real and 
positive within a selected spatial extent and normalized to 
a selected value to provide an updated estimate of h, then, 
using the updated estimate of h repeating the procedure 
until a selected criterion is met to provide a final 
estimate f as an output signal which is an improved 
estimate of f . 

13. The method of Claim 12 including the step of 
displaying the estimated signal f as an image on a display 
screen. 

14. The method of Claim 12 including, after the 
step of storing the signal data g, the step of segmenting 
the signal data g into overlapping segments in the spatial 
domain, and wherein the step of estimating comprises 
estimating for each segment of signal data g a selected 
bandlimit of the significant frequencies in the data, and 
including the additional steps of proceeding to the next 
segment of data and repeating the steps of the iterative 
procedure in the signal processor for each segment of data 
until estimated values of f and n are obtained for each 
segment of data, then recombining the estimates f for each 
segment of data to provide a combined estimated signal. 

15. The method of Claim 14 including the 
additional step of displaying the combined estimated signal 
as an image on a display screen. 

16. The method of Claim 12 wherein the initial 
estimate of the impulse response function h of the image 
acquisition system is a unit Dirac delta impulse function. 

17. The method of Claim 12 wherein the 
instrument is a confocal scanning microscope and wherein 
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1 the step of using the instrument to acquire a signal 

2 comprises scanning a specimen at a plurality of image 

3 planes at which the microscope is focused to develop a 

4 three dimensional array of image signal data in the spatial 

5 domain. 

1 18. The method of Claim 12 wherein the step of 

2 using the constrained estimate of h to obtain an estimate 

3 of f is carried out in accordance with the following steps 

4 at the kth iteration where and F k are the initial values 

5 or are determined at the k-1 iteration of n iterations: 

6 determining an estimator F k in the frequency 

7 domain as 

v XH k "(G-H k F k ) 
F k = F k + 

1 + XH k # H k 

8 for k = 0, 1, . . . . , n 

9 F Q = G, 

10 H p is the Fourier transform of a selected initial 

11 response function with selected bandlimits applied 

12 thereto; 

13 H k * is the complex conjugate of H k ; 

14 X is a selected parameter, and 

15 determining an estimate f k+l for the actual input 

16 signal as the inverse Fourier transform of F k for all 

17 values thereof which are real and positive within a 

18 selected spatial extent and as zero otherwise. 

1 19. The method of Claim 18 wherein the step of 

2 using the constrained estimate of f to obtain an estimate 

3 of h is carried out in accordance with the following steps 

4 wherein f k+1 is the estimate of f: 

5 determining the Fourier transform F k+1 of f k+1 ; 

6 determining an estimator H k as 
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where F" k+1 is the complex conjugate of F k+1 ; 



and 



applying the selected frequency bandlimits to H k ; 

determining an estimate h^, of the acquisition 
system instrument response function as the inverse Fourier 
transform of the bandlimited H k for all values thereof 
which are real and positive within a selected spatial 
extent and as zero otherwise and normalized to a selected 
value . 

20. The method of Claim 19 wherein the parameter 
X is selected based on the noise variance of the signal g, 
with lower values of X selected where high noise variance' 
is present and higher values of X selected where lower 
noise variance is present. 

21. The method of Claim 19 wherein the parameter 
X is selected to have a value in the range of 100 to 
20,000. 

22. The method of Claim 19 wherein the parameter 
X is selected to have a value of approximately 1,000. 

23. The method of Claim 12 wherein the variable 
of which the output signal g is a function is three 
dimensional . 

24. The method of Claim 22 wherein the 
instrument is a confocal microscope system and wherein the 
output signal g comprises digital data corresponding to 
light intensity values in a three dimensional array 
corresponding to light received from a three dimensional 
sample examined by the confocal microscope system. 
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1 25. An improved signal acquisition system having 

2 an instrument which provides an output signal g which is 

3 dependent on a variable of one or more dimensions, and a 

4 signal processor with a memory receiving the signal g in 

5 digital data form, the improvement comprising: 

6 (a) means in the signal processor for 

7 determining the transform G in the frequency domain of the 

8 output signal g; 

9 (b) means in the signal processor for 

10 receiving and storing in the memory selected frequency 

11 bandlimits in each dimension which are based on the range 

12 of signal frequency content exhibited by G in each 

13 dimension, and for receiving and storing an initial 

14 frequency domain estimate H 0 for the instrument having the 

15 selected bandlimits applied thereto; 

16 (c) means in the signal processor for 

17 iteratively determining estimates of the ideal input signal 

18 f to the instrument and of the actual instrument response 

19 function h, wherein the number of the last iteration is n, 

20 including at the kth iteration: 

21 means for determining an estimator F k in the 

22 frequency domain as 

F k s F k * DH k * (G-H k F k ) 

23 for k = 0 # 1, . • . . , n 

24 F Q = G, 

25 H 0 is determined as set forth above, 

26 H/ is the complex conjugate of H k ; 

27 D is a bounded linear operator; 

28 means for determining an estimate f k+1 for the 

29 actual input signal as the inverse transform to the spatial 
3 0 domain of F k f or all values thereof which are real and 

31 positive within a selected spatial extent and as zero 

32 otherwise; 

33 means for determining the transform to the 

34 frequency domain F k+ , of f k +ir 
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means for determining an estimator H k as 
H^H^DF/^G-H^) 

where F\ +1 is the complex conjugate of F k+l and D is a 
bounded linear operator; 

means for applying the selected frequency 
bandlimits to H k ; 

means for determining an estimate h^ of the 
acquisition system instrument response function as the 
inverse transform of the bandlimited H k for all values 
thereof which are real and positive within a selected 
spatial extent and as zero otherwise and normalized to a 
selected value; 

and wherein the means for iteratively determining 
continues to iterate until a selected criterion is met and 
for providing the final estimate as an output signal f 
which is an improved estimate of the ideal input signal f . 

26. The improved signal acquisition system of 
Claim 25 wherein 

1+\H/H k 



where X is a selected parameter, and 

l+XF kM F kM 



where X is a selected parameter • 

27. The improved signal acquisition system of 
Claim 25 wherein the initial estimate H 0 of the acquisition 
system instrument response function is the bandlimited 
frequency domain transform of the Dirac delta function. 
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1 28. The improved signal acquisition system of 

2 Claim 26 wherein the parameter X has a value in the range 

3 of 100 to 20,000. 

1 29. The improved signal acquisition system of 

2 Claim 25 wherein the parameter X is selected to have a 

3 value of approximately 1,000. 

1 30. The improved signal acquisition system of 

2 Claim 25 wherein the output signal g of the instrument is a 

3 function of a variable which is at least two dimensional, 

4 wherein the signal processor further includes means for 

5 storing data corresponding to the signal data g, for 

6 subdividing the signal data g into a plurality of sets of 

7 signal data, the signal data composing each set partially 

8 overlapping the signal data for each adjacent set, and for 

9 iteratively determining estimates of the ideal input signal 

10 f for each signal data set to determine a response 

11 function fi and estimated signal f for each signal data set, 

12 and for combining the estimated signal f for each signal 

13 data set and providing the combined estimated signals as 

14 the output signal. 

1 31. The improved signal acquisition system of 

2 Claim 2 5 wherein the variable of which the output signal g 

3 is a function is three dimensional. 

1 32. The improved signal acquisition system of 

2 Claim 25 wherein the acquisition system instrument is a 

3 confocal microscope and wherein the output signal g 

4 comprises digital data corresponding to light intensity 

5 values in a three dimensional array corresponding to light 

6 received from a three dimensional sample examined by the 

7 confocal microscope. 

1 33. The improved signal acquisition system of 

2 Claim 32 including a video screen and means for displaying 
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on the video screen to an operator a visual image of the 
output signal $ for data points corresponding to a single 
layer in a three dimensional array of data from the sample. 

34. The improved signal acquisition system of 
Claim 25 wherein the response function h^ is normalized to 
have a sum of values equal to one. 

35. The improved signal acquisition system of 
Claim 25 wherein the criterion for the last iteration is 
that a selected number of iterations is completed which is 
less than or equal to 10. 

36. The improved acquisition system of Claim 25 
wherein the transforms to the frequency domain are Fourier 
transforms and the transforms to the spatial domain are 
inverse Fourier transforms. 

37. An improved confocal microscope system 
having an optical microscope for examining a specimen, a 
laser scanning confocal microscope system coupled to the 
optical microscope to obtain optical image information 
therefrom and providing an output signal g representative 
of the microscope image, a signal processor connected to 
receive the output signal from the laser scanning confocal 
microscope system and to provide an output signal 
representing an optical image, and a display device 
connected to the signal processor to receive the output 
signal to display the image corresponding to the output 
signal from the signal processor, the improvement 
comprising: 

(a) means in the signal processor for 
determining the transform G in the frequency domain of the 
output signal g; 

(b) means in the signal processor for 
receiving and storing in the memory selected frequency 
bandlimits in each dimension which are based on the range 
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1 of signal frequency content exhibited by G in each 

2 dimension, and for receiving and storing an initial 

3 frequency domain estimate H Q for the conf ocal microscope 

4 system having the selected bandlimits applied thereto; 

5 (c) means in the signal processor for 

6 iteratively determining estimates of the ideal input signal 

7 f to the instrument and of the actual confocal microscope 

8 system response function h, wherein the number of the last 

9 iteration is n, including at the kth iteration: 

10 means for determining an estimator F k in the 

11 frequency domain as 

F k «F k + DH k # (G-H k F k ) 

12 for k « 0, 1, . . . n 

13 F 0 = G, 

14 1^ is determined as set forth above, 

15 H k * is the complex conjugate of H k ; 

16 D is a bounded linear operator; 

17 means for determining an estimate f k+1 for the 

18 actual input signal as the inverse transform to the spatial 

19 domain of F k for all values thereof which are real and 

20 positive within a selected spatial extent and as zero 

21 otherwise; 

22 means for determining the transform to the 

23 frequency domain F k+1 of t M ; 

24 means for determining an estimator H k as 

H^H^DF^G-H^,) 

25 where F* k+! is the complex conjugate of F k+I and D is a 

26 bounded linear operator; 

27 means for applying the selected frequency 

28 bandlimits to H k ; 

29 means for determining an estimate h t+l of the 

30 confocal microscope response function as the inverse 

31 transform of the bandlimited H k for all values thereof 

32 which are real and positive within a selected spatial 



WO 96/10794 



PCT/US95/12259 



- 65 - 

l extent and as zero otherwise and normalized to a selected 

> value; 

I and wherein the means for iteratively determining 

l continues to iterate until a selected criterion is met and 

> for providing the final estimate f a as the output signal to 

> the display. 

L 38* The confocal microscope system of Claim 37 

I wherein 

dh; - _ i 

1+XH/H k 

I where X is a selected parameter, and 



l+XP M r k#1 



L 39. The confocal microscope system of Claim 37 

> wherein the initial estimate H c of the confocal microscope 
l response function is the bandlimited frequency domain 

I transform of the Dirac delta function. 

L 40. The confocal microscope system of Claim 38 

I wherein the parameter X has a value in the range of 100 to 

I 20,000. 

41. The confocal microscope system of Claim 38 
J wherein the parameter X is selected to have a value of 

t approximately 1,000. 

42. The confocal microscope system of Claim 37 
! wherein the output signal g is a function of a variable 

\ which is at least two dimensional, wherein the signal 

l processor further includes means for storing data 

> corresponding to the signal data g, for subdividing the 
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1 signal data g into a plurality of sets of signal data, the 

2 signal data composing each set partially overlapping the 

3 signal data for each adjacent set, and for iteratively 

4 determining estimates of the ideal input signal f for each 

5 signal data set to determine a response function h and 

6 estimated signal function £ for each signal data set, and 

7 for combining the estimated signal f for each signal data 

8 set and providing the combined estimated signals as the 

9 output signal to the display device* 

1 43. The confocal microscope system of Claim 37 

2 wherein the variable of which the output signal g is a 

3 function is three dimensional. 

1 44. The confocal microscope system of Claim 37 

2 wherein the display device includes a video screen, and 

3 displays on the video screen to an operator a visual image 

4 of the estimated input function f for data points 

5 corresponding to a single layer in a three dimensional 

6 array of data from the sample. 

1 45. The confocal microscope system of Claim 37 

2 wherein the response function h k+l is normalized to have a 

3 sum of values equal to one. 

1 46. The confocal microscope system of Claim 37 

2 wherein the criterion for the last iteration is that a 

3 selected number of iterations is completed which is less 

4 than or equal to 10. 

1 47. The improved acquisition system of Claim 37 

2 wherein the transforms to the frequency domain are Fourier 

3 transforms and the transforms to the spatial domain are 

4 inverse Fourier transforms. 
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